Properties of the Fourier Transform¶

These demonstrations go beyond the content of the course. I don't expect you to understand what is going on in these snippets of code. But I show them to you so you see what sorts of things can be implemented by manipulating the Fourier coefficients.

Preliminaries¶

In [113]:
import numpy as np
from numpy.fft import fft, ifft, fft2, ifft2, fftshift, ifftshift
import matplotlib.pyplot as plt
%matplotlib inline
plt.ion()
from scipy import ndimage, misc
from scipy.signal import gaussian
In [114]:
from IPython.display import display, HTML
display(HTML("<style>.container { width:98% !important; }</style>"))

np.set_printoptions(edgeitems=30, linewidth=250, 
    formatter=dict(float=lambda x: "% .4f" % x))
In [115]:
def Gaussian2D(dims, std):
    '''
        Creates an array containing a 2-D Gaussian function, with mean
        at mu, and standard deviation std.

        Inputs:
          dims is a 2-vector containing the number of rows and columns
          std is the standard deviation
    '''
    nr, nc = dims
    rr, cc = np.mgrid[0:nr,0:nc]
    ctr = GetCentre(rr)
    temp1 = rr - ctr[0]
    temp2 = cc - ctr[1]
    temp1 = np.exp(-temp1**2/std)
    temp2 = np.exp(-temp2**2/std)
    blah = temp1*temp2
    return blah/sum(blah.flatten())

def GetCentre(f):
    return np.array(np.floor(np.array(np.shape(f))/2.), dtype=int)

Conjugate Symmetry¶

In [116]:
# If f is real-valued...
f = np.random.random(8)
print(f)

[ 0.6558  0.3726  0.1559  0.8886  0.5927  0.6509  0.9460  0.7658]
In [117]:
# Fourier coefs will be cojugate symmetric
F = np.fft.fft(f)
print(F)
print(F[2], F[-2])

[ 5.02823474+0.j         -0.22050572+0.90007346j  0.1466144 +0.63084762j  0.34663392-0.68017999j -0.32749292+0.j          0.34663392+0.68017999j  0.1466144 -0.63084762j -0.22050572-0.90007346j]
(0.146614397571887+0.6308476194085022j) (0.146614397571887-0.6308476194085022j)

Convolution¶

In [118]:
rr, cc = np.mgrid[-128:128, -128:128]

Try a circular kernel¶

In [119]:
%matplotlib inline
In [120]:
d = np.sqrt(rr**2 + cc**2)
circle = np.zeros([256,256])
circle_mask = d<15
circle[circle_mask] = 1.
circle = ifftshift(circle)
circle[0:22,0:5] = 1.
plt.figure(figsize=(8,4))
plt.subplot(1,2,1); plt.imshow(fftshift(circle), cmap='gray'); plt.axis('off'); plt.title('g, Shifted');
plt.subplot(1,2,2); plt.imshow(circle, cmap='gray'); plt.axis('off'); plt.title('g, Unshifted');
In [121]:
# Create an image with just a few non-zero pixels.
f = np.zeros([256,256])
f[50,100] = 1.0
f[48,110] = 2.
plt.figure(figsize=(8,8))
plt.imshow(f, cmap='gray'); plt.axis('off'); plt.title('f');
In [122]:
g = ifft2( fft2(f) * fft2(circle))

plt.figure(figsize=(8,8))
plt.imshow(abs(g), cmap='gray'); plt.axis('off');
plt.title(r'$(f \ast g)$');
In [123]:
f = plt.imread('pd.jpg')
f = f[:,:,0]
plt.imshow(f, cmap='gray'); plt.axis('off'); plt.title('An image f');
In [124]:
g = ifft2( fft2(f) * fft2(circle))

#plt.figure(figsize=(8,8))
plt.imshow(abs(g), cmap='gray'); plt.axis('off');
plt.title('g convolved with the image');

Or, try a Gaussian edge-detector kernel¶

In [125]:
# Make and edge detector
g = Gaussian2D([256,256], 20)
G = ifftshift( fft2(fftshift(g)) )
ramp = np.outer(range(256), range(256)) - 128**2
H = G * ramp*1.j
circle = np.real( ifft2( fftshift(H) ) )
plt.figure(figsize=[8,5]);
plt.subplot(1,2,1); plt.imshow(fftshift(circle), cmap='gray'); plt.axis('off'); plt.title('Edge Kernel (Shifted)');
plt.subplot(1,2,2); plt.imshow(circle, cmap='gray'); plt.axis('off'); plt.title('Edge Kernel (Unshifted)');
In [126]:
f = plt.imread('pd.jpg')
f = f[:,:,0]
In [127]:
F = fftshift( fft2( f ) )
In [128]:
plt.figure(figsize=[10,20])
plt.subplot(1,2,1)
plt.imshow(np.log(abs(F)+1), cmap='gray'); plt.axis('off'); plt.title('Shifted DFT (log-modulus)');
plt.subplot(1,2,2)
plt.imshow(f, cmap='gray'); plt.axis('off');
In [129]:
g = fftshift( ifft2( fft2(ifftshift(f)) * fft2(circle)) )

plt.figure(figsize=(16,8))
plt.subplot(1,2,1);
G = fftshift(fft2(ifftshift(g)))
plt.imshow(np.abs(G), cmap='gray'); plt.axis('off')
plt.title('DFT(Edge Kernel) x DFT(Image)');
plt.subplot(1,2,2); plt.imshow(abs(g), cmap='gray'); plt.axis('off');
plt.title('Edge Kernel convolved with image');

Filtering in the frequency domain¶

In [130]:
thresh = 10
rows, cols = np.shape(f)
ctr = np.floor(np.array(np.shape(f))/2) +1
rr, cc = np.mgrid[-ctr[0]:(rows-ctr[0]),
                  -ctr[1]:(cols-ctr[1])]
In [131]:
def PlotFiltImg(F, filt):
    F_filt = (filt)*F
    f_filt = ifft2( fftshift(F_filt) )
    plt.figure(figsize=[15,15])
    plt.subplot(2,2,1); plt.imshow(f, cmap='gray'); plt.axis('off');
    plt.subplot(2,2,3); plt.imshow(np.log(abs(F)+1), cmap='gray'); plt.axis('off');
    plt.subplot(2,2,4); plt.imshow(np.log(abs(F_filt)+1), cmap='gray'); plt.axis('off');
    plt.subplot(2,2,2); plt.imshow(np.real(f_filt), cmap='gray'); plt.axis('off');
In [132]:
filt_mask = abs(cc)<thresh
filt = np.zeros(np.shape(F))
filt[filt_mask] = 1.0
PlotFiltImg(F, filt)
In [133]:
filt_mask = abs(rr)<thresh
filt = np.zeros(np.shape(F))
filt[filt_mask] = 1.0
PlotFiltImg(F, filt)
In [134]:
filt_mask = np.sqrt(rr**2 + cc**2)<thresh
filt = np.zeros(np.shape(F))
filt[filt_mask] = 1.0
PlotFiltImg(F, filt)

Applying that filter in the frequency domain is the same as convolving the image with the kernel below, attained by taking the IDFT of the filter kernel.

In [135]:
Filt = fftshift(ifft2(ifftshift(filt)))
plt.figure()
plt.subplot(1,2,1); plt.imshow(filt, cmap='gray'); plt.title('DFT of filter'); plt.axis('off')
plt.subplot(1,2,2); plt.imshow(np.log(abs(Filt)+1), cmap='gray'); plt.axis('off'); plt.title('Spatial filter');

Scaling¶

You can scale an image using the DFT by padding in cropping in the two domains (time/space, and frequency).

In [136]:
# We will magnify an image by padding its DFT coefficients, and then cropping the image
# that results from the IDFT of those padded coefs.
scale_factor = 1.5
padamt = int(rows*(scale_factor-1.)/2.)
G = np.pad(fftshift(fft2(f)), (padamt,), mode='constant')
g512 = ifft2( ifftshift(G) )
if padamt!=0:
    g = g512[padamt:-padamt,padamt:-padamt]
else:
    g = g512
In [137]:
plt.figure(1,figsize=[15,15])
plt.clf()
plt.subplot(1,2,1)
plt.imshow(f, cmap='gray'); plt.title('Original')
plt.subplot(1,2,2)
plt.imshow(np.real(g), cmap='gray'); plt.title('Scaled x'+str(scale_factor));

Fourier Shift Theorem¶

You can shift an image by adding a linear 'ramp' to the phase portion of its DFT coefficients.

In [139]:
dr = 20
dc = 12
ramp = -2.*np.pi*1j*(cc*dc + rr*dr)/rows
ramp = np.roll(np.roll(ramp,-1, axis=0),-1, axis=1)
wave = np.exp(ramp)
plt.figure()
plt.imshow(np.imag(ramp))
plt.title('Ramp in imaginary part');
In [140]:
# Make sure the centre pixel of the ramp has a phase of zero.
ctr = int(np.floor(256/2))
print('This should be zero -> '+str(ramp[ctr,ctr]))
print('This is the DC of our image -> '+str(F[ctr,ctr]))  # This should be the DC of our image.

This should be zero -> -0j
This is the DC of our image -> (2491737+0j)
In [141]:
G = wave * fftshift(fft2(f))
g = ifft2(ifftshift(G))

plt.figure(1,figsize=[15,15])
plt.subplot(1,2,1)
plt.imshow(f, cmap='gray'); plt.title('Original')
plt.subplot(1,2,2)
plt.imshow(np.real(g), cmap='gray');
plt.title('Shifted (down,right)=('+str(dr)+','+str(dc)+')');

Pendulum synchronization video (YouTube)

Phase Correlation¶

You can detect a shift between two similar images by looking for a ramp in their phase difference. Or, you can take that phase difference and take the IDFT of it, and look for a spike. The offset of the spike from the origin gives the relative shift.

In [143]:
plt.figure(figsize=(8,8));
g = ndimage.shift(f,[dr,dc]) + np.random.normal(0.0, 5, size=np.shape(f))
G = fft2(g)
F = fft2(f)
H = G/F
plt.imshow(np.real(fftshift(H)), cmap='gray', vmin=-4, vmax=4);
In [144]:
%matplotlib notebook
plt.figure(figsize=(8,8));
h = ifft2(H)
plt.imshow(np.real(h), cmap='gray');
In [145]:
# Where does that lonely pixel occur?
coords = np.unravel_index(np.argmax(np.real(h)), h.shape)
print('Bright pixel at '+str(coords))
print('True shift was '+str((dr,dc)))

Bright pixel at (20, 12)
True shift was (20, 12)

Rotation¶

You can also rotate an image by simply applying the same rotation to the frequency domain.

In [146]:
theta = 30
f = plt.imread('pd.jpg')
f = np.array(f[:,:,0], dtype=float)
In [147]:
padamt = int(rows/2)
#padamt = 0
In [148]:
# Padding the image first makes the frequency domain smoother
ff = f
ff = np.pad(f, (padamt,), mode='constant')
FF = fftshift( fft2( ifftshift(ff) ) )
In [149]:
# Apply rotation to the Fourier coefficients
Gr = ndimage.rotate(np.real(FF), theta, reshape=False, order=1)
Gi = ndimage.rotate(np.imag(FF), theta, reshape=False, order=1)
G = Gr + 1j*Gi
In [150]:
g = fftshift( ifft2( ifftshift(G) ) )
if padamt!=0:
    g = g[padamt:-padamt,padamt:-padamt]
In [152]:
%matplotlib inline
plt.figure(figsize=[15,15])
plt.subplot(1,2,1)
plt.imshow(f, cmap='gray'); plt.title('Original'); plt.axis('off')
plt.subplot(1,2,2)
plt.imshow(np.real(g), cmap='gray'); 
plt.title('Rotated by '+str(theta)+' degrees'); plt.axis('off');

Fourier Projection Theorem¶

The Fourier projection theorem is the basis for reconstruction in Computed Tomography (CT) scans.

In [153]:
f = plt.imread('ct.jpg')
f1 = np.array(f[:,:,0], dtype=float)
In [154]:
plt.figure(1, figsize=(10,10))
plt.clf()
plt.subplot(2,2,1)
plt.imshow(f1, cmap=plt.cm.gray);

px = np.sum(f1,axis=0)
plt.subplot(2,2,3)
plt.plot(px)
plt.axis('tight');

F = fftshift(fft2(f1))
plt.subplot(2,2,2)
Fnice = np.log(abs(F)+1)
plt.imshow(Fnice, cmap=plt.cm.gray);

ctr = GetCentre(F)
Fx = F[ctr[0],:]
plt.plot([1,254],[128,128], 'r--')
plt.axis([0,255,255,0])

plt.subplot(2,2,4)
plt.plot(np.log(abs(Fx)+1), 'r')
plt.axis('tight');

Px = fftshift( np.fft.fft(px) )
#plt.subplot(2,2,4)
plt.plot(np.log(abs(Px)+1), 'b--');

Aliasing in images¶

In [155]:
x, y = np.mgrid[0:256,0:256]
f = np.sin(2.*np.pi*x/256*40)
#f = ndimage.imread('pd.jpg')
plt.figure(figsize=(15,7))
plt.clf()
plt.subplot(1,3,1)
plt.imshow(f, cmap='gray');
c = plt.axis()
plt.title('Original');

g = f[::4,::4]  # every fifth sample
plt.subplot(1,3,2)
plt.imshow(g, cmap='gray'); #plt.axis(c);
plt.title('Subsampled');

plt.subplot(1,3,3);
plt.imshow(g, cmap='gray'); plt.axis(c);
plt.title('Subsampled');
In [156]:
q = plt.imread('bricks.jpg')
q = np.array(q[:,:,0], dtype=float)
nr, nc = np.shape(q)
In [157]:
plt.figure(figsize=(15,7));
plt.subplot(1,3,1)
plt.imshow(q, cmap='gray', interpolation='none')
plt.title('Original'); c = plt.axis();

q_small = q[::5,::5]
qsr, qsc = np.shape(q_small)
plt.subplot(1,3,2)
plt.imshow(q_small, cmap='gray', interpolation='nearest')
plt.axis(c);
plt.title('Subsampled, with aliasing');

plt.subplot(1,3,3)
plt.imshow(q_small, cmap='gray', interpolation='nearest')
plt.title('Subsampled, with aliasing');
In [158]:
plt.figure(figsize=(15,7));
plt.subplot(1,3,1)
plt.imshow(q, cmap='gray', interpolation='none')
plt.title('Original');

q_small = q[::5,::5]
qsr, qsc = np.shape(q_small)

g = Gaussian2D(q.shape, 10)   # Create filter
Q = fftshift( fft2(q) )       # DFT of the filter

G = fftshift( fft2(fftshift(g)) )
h = ifft2(ifftshift(Q*G))     # Apply filter in the freq domain

plt.subplot(1,3,2)
plt.imshow(np.real(h[::5,::5]), cmap='gray', interpolation='nearest')
plt.axis(c);
plt.title('Smoothed and subsampled (no aliasing)');

plt.subplot(1,3,3)
plt.imshow(np.real(h[::5,::5]), cmap='gray', interpolation='nearest')
plt.title('Smoothed and subsampled (no aliasing)');
In [ ]: